{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b8674439",
   "metadata": {},
   "source": [
    "# Reactive Transport\n",
    "\n",
    "> This notebook is still a work in progress\n",
    "\n",
    "Incorporating chemical reactions or source terms in general is an important feature. This notebook will explain how OpenPNM implements this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "979d4ed7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openpnm as op\n",
    "import numpy as np\n",
    "op.visualization.set_mpl_style()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "762000f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "net : <openpnm.network.Cubic at 0x253943b3860>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.coords                                                         9 / 9\n",
      "  3  throat.conns                                                      12 / 12\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.surface                                                            8\n",
      "  3  throat.surface                                                          8\n",
      "  4  pore.left                                                               3\n",
      "  5  pore.right                                                              3\n",
      "  6  pore.front                                                              3\n",
      "  7  pore.back                                                               3\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "pn = op.network.Cubic(shape=[3, 3, 1])\n",
    "print(pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ced2e267",
   "metadata": {},
   "source": [
    "And let's also create a ``Phase`` object to store the diffusive conductance values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "e929f949",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph = op.phase.Phase(network=pn)\n",
    "ph['throat.diffusive_conductance'] = np.random.rand(pn.Nt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4794959b",
   "metadata": {},
   "source": [
    "## Using the source term functionality in OpenPNM\n",
    "\n",
    "First, we'll illustrate the process of adding a source term to a diffusion simulation, then we'll look behind the scenes at how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0763cfd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd = op.algorithms.FickianDiffusion(network=pn, phase=ph)\n",
    "fd.set_value_BC(pores=pn.pores('left'), values=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8a665fe",
   "metadata": {},
   "source": [
    "OpenPNM's `models` module contains a variety to pre-written source terms. Let's use a standar reaction, also known as a power-law:\n",
    "\n",
    "$$ y = A_0 X^{A_1} + A_3 \\rightarrow r_A = kC_A^b$$\n",
    "\n",
    "where $C_A$ is the concentration of the diffusion species $A$, $k$ and $b$ are the kinetic constant (actually a strong function of temperature) and $b$ is the reaction order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bedb6b58",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = op.models.physics.source_terms.power_law"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed4fce74",
   "metadata": {},
   "source": [
    "This model is added to the `ph` object, since it is a function of the phase properties, specifically concentration and temperature if the reaction is not isothermal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "713445af",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph.add_model(propname='pore.reaction',\n",
    "             model=f,\n",
    "             X='pore.concentration',\n",
    "             A1=1e-5,\n",
    "             A2=1,\n",
    "             A3=0.0,\n",
    "             regen_mode='deferred')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0d838ff",
   "metadata": {},
   "source": [
    "Let's unpack the above code block:\n",
    "- The model is being stored under `'pore.reaction'`, which means that the `FickianDiffusion` algorithm will call this model to retrieve the reaction rate at the given conditions. However, we still have to tell `fd` about this reaction.\n",
    "- The model uses the values of `'pore.concentration'` to compute $r_A$, so we specify to use that array for `X`.\n",
    "- Since `'pore.concentration'` does not exist yet this model will throw an error if run, so we set the `regen_mode` to `'deferred'`. \n",
    "- `A1`, `A2` and `A3` can be either scalars or dictionary keys to numpy arrays on the `ph` object. This is possible since all OpenPNM dict object return any non-string keys directly back, so `fd[1.0]` returns `1.0`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22fe6678",
   "metadata": {},
   "source": [
    "Now we add this model to our algorithm, to tell the algorithm about it and also in which pores it applies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ec2ed8e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd.set_source(propname='pore.reaction', pores=pn.pores('right'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d88299",
   "metadata": {},
   "source": [
    "The act of adding a source term creates a new array on the algorithm as can be seen below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "56eff1c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "fick_02 : <openpnm.algorithms.FickianDiffusion at 0x253953195e0>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.bc.rate                                                        0 / 9\n",
      "  3  pore.bc.value                                                       3 / 9\n",
      "  4  pore.concentration                                                  9 / 9\n",
      "  5  pore.initial_guess                                                  9 / 9\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.all                                                                9\n",
      "  3  throat.all                                                             12\n",
      "  4  pore.source.reaction                                                    3\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(fd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399c5437",
   "metadata": {},
   "source": [
    "The `'pore.source.reaction'` array contains a `boolean` array where `True` values indicate the locations where the source term should be applied.  Had we added a second reaction, called `'pore.another_reaction'`, the above printout would have an entry for `'pore.source.another_reaction'`.  This is designed to make it easy to get a list of all source terms and their location using `fd['pore.source']` which would return a `dict` with `'reaction'` and `'another_reaction'` as keys, and two `boolean` arrays as values indicating where each reaction was applied."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65cc2581",
   "metadata": {},
   "source": [
    "So now the algorithm is ready to be run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "35d876a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                           \r"
     ]
    }
   ],
   "source": [
    "fd.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0126b8dc",
   "metadata": {},
   "source": [
    "Let's look at which arrays are created on both `fd` and `ph`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "1b70d986",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "fick_02 : <openpnm.algorithms.FickianDiffusion at 0x253953195e0>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.bc.rate                                                        0 / 9\n",
      "  3  pore.bc.value                                                       3 / 9\n",
      "  4  pore.concentration                                                  9 / 9\n",
      "  5  pore.initial_guess                                                  9 / 9\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.all                                                                9\n",
      "  3  throat.all                                                             12\n",
      "  4  pore.source.reaction                                                    3\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(fd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4214e2c",
   "metadata": {},
   "source": [
    "Above we can see that `'pore.concentration'` was created, which contains the solution of the transport problem. `'pore.initial_guess'` was also created since the algorithm was aware that a source term was added, so iteration of the problem would be required hence an initial guess was needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f51dbc6a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "phase_02 : <openpnm.phase.Phase at 0x253943cdc20>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.temperature                                                    9 / 9\n",
      "  3  pore.pressure                                                       9 / 9\n",
      "  4  throat.diffusive_conductance                                      12 / 12\n",
      "  5  pore.concentration                                                  9 / 9\n",
      "  6  pore.reaction.S1                                                    9 / 9\n",
      "  7  pore.reaction.S2                                                    9 / 9\n",
      "  8  pore.reaction.rate                                                  9 / 9\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  2  pore.all                                                                9\n",
      "  3  throat.all                                                             12\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(ph)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b7d117",
   "metadata": {},
   "source": [
    "We can see that `'pore.concentration'` was *also added* to this object. This was necessary since the source-term models needed access to this array in order to compute the reaction rate. So the `fd` algorithm writes this array to `fd` and *then* regenerates the source term model(s). Recall that we set `regen_mode='deferred'` when assigning this model, so the algorithm is the first to run this model. \n",
    "\n",
    "We can also see `'pore.reaction.S1'`, `'S2'` and `'rate'` which are all values returned by the model. We'll explain what `'S1'` and `'S2'` mean in the following section, and `'rate'` is self explanatory. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82450d01",
   "metadata": {},
   "source": [
    "## How the source-term functionality works\n",
    "\n",
    "OpenPNM uses the method described by Patankar [ref] to linearize the source-terms.  As will be shown below, this method computes the slope and intercept of non-linear term at the current value of the quantity being computed.  This is the origin of the `'S1'` (the slope) and `'S2'` (the intercept) arrays on the `ph` object above.\n",
    "\n",
    "Determination of `'S1'` and `'S2'` will be outlined below, but let's first illustrate how they are used. Consider the mass balance around pore 6, which connects to pore 3 and 7:\n",
    "\n",
    "$$ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = r_A $$\n",
    "\n",
    "Meaning the net diffusive flux into pore 6 from pores 3 and 7 is equal to the rate of consumption in pore 6.  We can further write the rate value as:\n",
    "\n",
    "$$ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = k C_{A,6}^2 $$\n",
    "\n",
    "Now we can see there will be a problem solving this equation since a non-linear concentration term has appeared.  *If* we could just write the reaction term as a linear expresssion we would be able to proceed:\n",
    "\n",
    "$$ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = S_1 C_{A,6} + S_2 $$\n",
    "\n",
    "We could then group all the linear terms as follows:\n",
    "\n",
    "$$ g_{3,6} C_{A,3} + (-g_{3,6} - g_{7,6} - S_1) C_{A,6} + g_{7,6} C_{A,7} =  S_2 $$\n",
    "\n",
    "The above equation (or rather the system of equations to which is belongs) is then solved using standard numerical solvers to find all the unknown values of $C_A$. Now we just need to find `'S1'` and `'S2'`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44172700",
   "metadata": {},
   "source": [
    "The process starts with the standard definition of Newton's method:\n",
    "\n",
    "$$ S'=\\frac{(S-S^*)}{(x_i-x_i^* )} $$\n",
    "\n",
    "This can be rearranged to solve for $S$:\n",
    "\n",
    "$$ S=S^*+(S')^*⋅(x_i-x_i^* ) $$\n",
    "\n",
    "Where $S'^*$ is is the derivative of the function at $x_i^*$. Since $x_i^*$ is known (as either an initial guess or a updated value from each iteration), we can treat it and $S'^*$ as constants and lump them together as follows:\n",
    "\n",
    "$$ S= (S' )^* x_i+(S^*-(S' )^* x_i^* )=S_1 x_i+S_2 $$\n",
    "\n",
    "So both $S$ values are computed from known values, and the result is a linear equation that is a function of the unknown $x_i$, which can be entered into the system of equations as mentioned in the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "156d1cb3",
   "metadata": {},
   "source": [
    "Now that we have seen how `'S1'` and `'S2'` are defined and computed, we can  revisit how this is implemented in OpenPNM. We the `run` method is called on the `fd`, several hidden methods are called, these include:\n",
    "\n",
    "- `_build_A`\n",
    "- `_build_b`\n",
    "- `_apply_bcs`\n",
    "- `_apply_sources`\n",
    "\n",
    "Let see how the coefficient and RHS matrix evolve as these are each called."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "f1f15bda",
   "metadata": {},
   "outputs": [],
   "source": [
    "fd = op.algorithms.FickianDiffusion(network=pn, phase=ph)\n",
    "fd.set_value_BC(pores=pn.pores('left'), values=1.0)\n",
    "fd.set_source(propname='pore.reaction', pores=pn.pores('right'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "39e1abe9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.49827553 -0.62279691  0.         -0.87547862  0.          0.\n",
      "   0.          0.          0.        ]\n",
      " [-0.62279691  1.49321788 -0.77524129  0.         -0.09517969  0.\n",
      "   0.          0.          0.        ]\n",
      " [ 0.         -0.77524129  1.22425256  0.          0.         -0.44901127\n",
      "   0.          0.          0.        ]\n",
      " [-0.87547862  0.          0.          1.43888372 -0.21067434  0.\n",
      "  -0.35273075  0.          0.        ]\n",
      " [ 0.         -0.09517969  0.         -0.21067434  0.92271523 -0.18960804\n",
      "   0.         -0.42725315  0.        ]\n",
      " [ 0.          0.         -0.44901127  0.         -0.18960804  1.06621712\n",
      "   0.          0.         -0.42759781]\n",
      " [ 0.          0.          0.         -0.35273075  0.          0.\n",
      "   0.40896966 -0.05623891  0.        ]\n",
      " [ 0.          0.          0.          0.         -0.42725315  0.\n",
      "  -0.05623891  0.5941874  -0.11069534]\n",
      " [ 0.          0.          0.          0.          0.         -0.42759781\n",
      "   0.         -0.11069534  0.53829315]]\n"
     ]
    }
   ],
   "source": [
    "print(fd.A.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b710638",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
